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WSLD operators: A class of fourth order difference approximations for 
space Riemann-Liouville derivative 
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Because of the nonlocal properties of fractional operators, higher order schemes play more important role 
in discretizing fractional derivatives than classical ones. The striking feature is that higher order schemes 
of fractional derivatives can keep the same computation cost with first-order schemes but greatly improve 
the accuracy. Nowadays, there are already two types of second order discretization schemes for space 
fractional derivatives: the first type is given and discussed in [Sousa & Li, arXiv: 1 109.2345 , Chen & 
Deng, arXiv: 1304.3788 , Chen et al., Appl. Numer. Math., 70, 22-41]; and the second type is a class of 
schemes presented in [Tian et al., arXiv:1201.5949|. The core object of this paper is to derive a class 
of fourth order approximations, called the weighted and shifted Lubich difference (WSLD) operators, 
for space fractional derivatives. Then we use the derived schemes to solve the space fractional diffusion 
equation with variable coefficients in one-dimensional and two-dimensional cases. And the unconditional 



stability and the convergence with the global truncation error 
numerically verified. 



(t +/; ) are theoretically proved and 



Keywords: Fractional diffusion equation; Weighted and shifted Lubich difference operators; Numerical 
stability; Convergence 



1. Introduction 



In rec ent decades, fractional operators have been playing more and more important roles llDiethelm 
(I2OI0I) ]. e.g., in mechanics (theory of viscoelasticity and viscoplasticity), (bio-)chemistry (modelling 
of polymers and proteins), electrical engineering (transmission of ultrasound waves), medicine (mod- 
elling of human tissue under mechanical loads), etc. Efficiently solving the fractional partial differential 
equations (PDEs) naturally becomes an urgent topic. Because of the nonlocal properties of fractional 
operators, obtaining the analytical solutions of the fractional PDEs is more challenge or sometimes even 
impossible; or the obtained analytical solutions are less valuable (expressed by transcendental functions 
or infinite series). Luckily, some important progr ess has been made for numeric a lly solving the frac - 
tional PDEs by finite difference methods, e.g., see OMeerschaert & TadieranI (120041) ; ISousa & Lil (1201 lb ; 
Sun & Wul (l2006l) ; lTian et al\ (I2012h : l^ist3(l2006l) : lzhuang et al\ (l2009h ]. 

In solving space fractional PDEs, high order finite difference schemes display more striking ben- 
efits because most of the time they can use the same computational cost with first order scheme but 
greatly improve the accuracy. For example, comparing with first order difference scheme which may 
have the matrix algebraic equation (/— A)m"+' — u" +b"^\ the high order scheme has the matrix al- 
gebraic equation (/— A)m"+' ~ {I + B)u" +b"^^i^. The three matrices A, A and B are all Toeplitz-like 
and have completely same structure, and the computational count for matrix vector multiplication is 
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^(Nlo ^N), then the com putational costs for solving the two matrix algebraic equations are almost the 
same liChen ef ad(l20l2l) 1. 

Nowadays, we notice that there exist two types of second order discretization schemes for space frac- 
tional derivatives. The idea of the first type is to combine the centered difference scheme of second clas- 
sical d erivative with piecewise linear polynomial approximation of the fractional integral. ISousa & Li 
( 2011) firstly use the idea to obtain the second order approximation in infinite domain. The paper 
OChen & Deng (1201 ll) 1 detailedly analyzes the effectiveness of the approximation in finite domain. And 
this discretization is also effectively used to solve the time-space Capuo-Riesz fractional diffusion equa- 
tion lichen et al.\ (l2013h l. The second type of second order approximation is in fact a class of second 
order discretization, which are obtained by assembling the Griinwald difference operators with differ- 
ent weights and shifts. This class of approxima tions are detailedl y discussed and successfully applied 
to solve space fractional diffusion equations in BTian et al.\ ( l2012h l. and called WSGD operators there. 
Both of the two types of the operators have completely same stru c ture, and the real parts of the eigen- 
values of the matrixes are less than 0, see ODeng & ChenI (120 13h : iTian et al\ (120121) 1. So they can be 
efficiently used to solve space fra ctional PDEs. 

Based on Lubich's operator [Lubich] (119861) 1. this paper derives a class of fourth order approxi- 
mations for space fractional derivatives, termed the weighted and shifted Lubich difference operators 
(WSLD operators). Using the fractional linear multistep methods, iLubichI (119861) obtains the L-th or- 
der (L ^ 6) approximations of the a-th derivative (a > 0) or integral (a < 0) by the corresponding 
coefficients of the generating functions 5"{Q, where 



5"(0=fEj(i-c)'j 



(1.1) 



For a — 1, the scheme reduces to the classical (L+ 1) -point backward difference formula OHenrici 
(1 19621) 1 . For L = 1 , the scheme (ILlt corresponds to the standard Griinwald discretization of a-th deriva- 



tive with first order accuracy; unfortunately, f or the time dependent equations the difference scheme is 
unstable. But lMeerschaert & TadjeranI (|2()04|) successfully cir cumvent this difficulties by the so-called 
shifted Griiwald formulae. Taking L = 2. ICuesta et al\ ( 120061) discuss the convolution quadrature time 
discretization of fractional diffusion-wave equations; when applying the discretization scheme to space 
fractional operator with a G (1,2) for time dependent problem, the obtained scheme is also unstable, 
since the eigenvalues of the matrix corresponding to the discretized operator are greater than one. If 
using the shifted Lubich's formula, it reduces to the first order accuracy (detailed description is given in 
Section 2). This paper weights and shifts Lubich's operator to obtain a class of fourth order discretiza- 
tion schemes, which are effective for time dependent problem. Then we use the fourth order schemes to 
solve the following two-dimensional fractional diffusion equation with variable coefficients. 



' du(x,y,t) 



'd+{x,y)^^Dy{x,y,t)+d-{x,y):,D"u{x,y,t) 



+ e+{x,y)yi^D^u{x,y,t)+e-{x,y)yD^^u{x,y,t)+f{x,y,t), 
u{x,y,0)^uo{x,y), for {x,y)eQ, 
, u{x,y,t) —0, for {x,y,t) e dQ x {0,T], 



(1.2) 



in the domain Q = (x^jXr) x {yL,yR), Q < t ^ T, where the orders of the fractional derivatives are 
1 < a,j3 < 2 and f{x,y,t) is a source term, and all the variable coefficients are nonnegative. The left 
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and right Riemann-Liouville fractional deriva t ives of the function m(x) on [jc^.j^a], —°°^xl<xr^°o 
are, respectively, defined by [iPodlubnvl(ll999l) : ISanik:o et a/.l ( ll993b l 



D'}u{x) 



xl'^x 



and 



^Dl<^) 



r{2-a)dx^ 



1 



r{x-^y-"u{^)d^, 

J Xl 



{^-x)'-%i^)d^. 



(1.3) 



(1.4) 



r{2-a)dx^ . 

The outline of this paper is as follows. In Section 2, we derive a class of fourth order approxi- 
mations for space fractional Riemann-Liouville derivatives, being effective in solving space fractional 
PDEs. In Section 3, the full discretization schemes of one-dimensional case of ( ll.2l i and ( 11.2b itself are 
presented. Section 4 does the detailed theoretical analyses for the stability and convergence of the given 
schemes. To show the effectiveness of the algorithm, we perform the numerical experiments to verify 
the theoretical results in Section 5. Finally, we conclude the paper with some remarks in the last section. 



2. Derivation of a class of fourth order discretizations for space fractional operators 

In the following, we derive a class of fourth order approximations for Riemann-Liouville fractional 
derivatives, and prove that they are effective in solving space fractional PDE, i.e., all the eigenvalues of 
the matrixes corresponding to the discretized operators have negative real parts. 



2.1 Derivation of the discretization scheme 

Taking L = 2, for all | ^| ^ 1, Eq. dl.lb can be recast as 



2 ^2^ 



with k = m + n, and 



rf = (-i)M^ L3 



3\" 1 

2) (i-c)"(i-3cr 



2) K-D" "c"-L 



«=0 



a 



k=0 



n=0 
a fk 



£(-!)" (!^ 



m=0 



a k 



a 

k — m 



a 



ym+n 



a k 



2^ '^ omSk—m^ 
m=0 



(2.1) 



(2.2) 



a 



where ^" = (— 1)*^ ( "T ) are the coefficients of the power series of the generating function (1 — Q", and 
they can be calculated by the following recursively formula 

a + l 



g^^l, rf= 1- 



rf-i, k^i- 



(2.3) 
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If a < 0, {?"}^=o correspond to the coeffic ients of the 2nd orde r convolution quadrature for the approx- 
imation of fractional integral operator [see. lCuesta et g/.l (120061) 1. 



Lemma 2.1 The coefficients in (12.21 ) with a € (1,2) satisfy the following properties 



9S=l^l >0; 

„_ /3\"a(a-l)(64a2-176a + 123) 
"^'-[ij 486 ''"' 

3\"2a(a-l)(2-a)(64a2- 208a + 183) 



?f 



ft" 



3\"4a 



<0; 

3^|"4a(a -l)(7-8a) 
^1 



27 3 
a 



<0; 



^^ V2y 3645 

Proof. Taking i^ = 1, it is easy to check that 



>0; 



k=0 



a fk 



k=Q k=Q 



2-2C+2C' 



n 

We first introduce two lemmas, which will be used to prove that the several classes of derived 
discretization schemes are 2nd, 3rd, and 4th order convergent, respectively. 

Lemma 2.2 (JErvin & RoopI (I2OO6I) ) Let a > 0, m G Cq{Q), Q cR, then 



^(-ooDXx)) = {-ico)''u{co) and ^(;,D^m(x)) = (/«)"«(«), 
where ^ denotes the Fourier transform operator and u{(o) = ^{u), i.e., 



u{(0) 



'u{x)dx. 



Lemma 2.3 Let m, _ooZ)"+'m(x) (or _ooZ)"+^m(x)) with a G (1,2) and their Fourier transforms belong 
to Li (R) when p j^Q {01 p = 0); and denote that 



i^"«W = 7^ E qfu{x-{k-p)h), 
" k=0 



(2.4) 



where q" is defined by ( I2.2l l and p an integer. Then 



and 



„DX;c) = i4pM(^) + ^(/i), P^O, 
.DXx) = lA"m(x) + &{h^), p = 0. 
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Proof. From (I2.2l i and k = m + n,we obtain 



k=0 

k=Q ^ ' 



a —iaph 



h-'^e 



{-icof 



E(-i) 



-m)''e"' 



-impli 
1- 



n=0 
1- 



" f I Joynh 



1 \ '" / \ 



,lfi)/! \ " 



-icoh 



;n=0 



l+^ll- 



iCDmh 



ii{co) 



iah 



u{co) 



l + -{l-e-')] uico), 



with z = —icoh. It is easy to check that 



1-, 






6^ -4'' 



24 



48 



and 



1 + 2(1-^"^)) = 



" ' a a(a-3) 2 a(a2-9a + i2)3 _^ .; 

1 + -z + ^ — -z^ + — -t + G{z) 



then we have 

'\-e 



l + 2(^-' 



48 



3/72 -2a 2 2/73 + a(3-4/7) 3 ^,4, ^^ ^^ 

l+pz+-^-^ z^ + -^ 77; —z'^e(z'). (2.5) 

O 12 



Therefore, from Lemma l2n we get 



^aA««)(«) = ^(_D««(x)) + 0(0)), 



where (co) = (-/«)« [pz + ^£^22 _^ V+a(3-4p) ^3 ^ ^(^4)^ ^(^)^ Tjjgjj tijgjg g^i^jg 



|0(a))|^?l/a)|"+'|M(a))|-/i, /?^0, 
|0(«)| s;c|^■«|"+V(«)|•/^^ P = 0- 



Hence 



|_oeDXx) ~ lA«m(x)| = |(/)(x)| s; ^ L l^(«)l^-« 



^(/i2), p = 0. 



D 
In the following, we present the approximation operators for Riemann-Liouville derivative and prove 
that they have 2nd, 3rd, and 4th order truncation errors. 
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Theorem 2. 1 (Second order approximations for left Riemann-Liouville derivative) Let m, _ooZ)"^^m (x) 
with a e (1,2) and their Fourier transforms belong to Li(M). Denote that 



2lA"^u{x) = WpLApU(x)+Wc,LA"u{x), 
where tA" ^A" are defined by ( I2.4l l. w^ = -^, w,, = -^, p ^ q, and p, ^ are integers. Then 

Proof. From the proof of Lemma |23] we have 



(2.6) 



^aA>)(a)) = (-/a))« 



and 



Then there exists 



l+pz- 



l+qz- 



3p2-2a 2 2p3 + a(3-4/?) 



12 



zV^(z4) «(«) 



3^2 -2a 2 2cy^ + a(3-4^) 
6 ^ ^ 12 



z^ + ff{z^) u{co). 



3pq + 2a 9 2pq{p + q) — 3a 

1 ; Z" -:: Z 



and by the similar way to the proof of Lemma |2J1 we get 



6 ^^- """ 12" ^' + ^(2')J«(®), 



n 



Theorem 2.2 (Third order approximations for left Riemann-Liouville derivative) Let u, _aoZ)" m(x) 
with a E (1,2) and their Fourier transforms belong to Li(]R). Denote that 



3L^"?.M«W = Wp^^2LA"^gU{x) +Wr,s2LAZAx), 



(2.7) 



_ 3r.v+2a ,„ _ 3pq+2a 



where 2lA".,, and 2lA". are defined by ^M, Wp^g = j(^^y w,-^ = ^^^, « 7^ /??, and p, q, r, i are 
integers. Then 

_^««(^) - 3lA,^,„,,«(x) + ^(/l^). 

Proof. By the proof of Theorem l2.1l we have 

3pq + 2a 2 2pq{p + q) — 3a 3 



and 



^(2LA»(a)) ==(-/«)« 



^(2lA« «)(«) = i-iCoY 



1 



12 



z^ + ^(z^) m(«) 



3rj' + 2a , 2rs(r + s) -3a 
1 : Z ^ — ^ i 



g -r -^^^ z^ + ^(z^) «(«). 



Then there exists 



^(3LA«,,,,,«)(a)) 



(0) 



6/7^r,s-(r + ,s--/:)-^)+4a[ri(r + 5)-/5^(p + ^)]+9a(ri-/?^) 3 , „. 4. ^. . 
' + 36(..,-p,) ^ +'''^ ' "<'»'• 
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and by the similar way to the proof of Lemma lZlJ we get 



D 



Theorem 2.3 (Fourth order approximations for left Riemann-Liouville derivative) Let u, -ooD^^^u{x) 
with a G (1,2) and their Fourier transforms belong to Li(E). Denote that 



where sl^? ^ ^^. and 3^A^-^^ are defined by ( |2.71 >; and 



''p,q,r,s 



^ p,q,r,s 



^p.qj\s ^,q,r,s ^p,q,r.s ^ p,q,}\s 

^"pJj.T.^^p.qj.s 
^"p^J^^p.q.r.s ^p,q,r,s^^^^ 



(2.8) 

(2.9) 
(2.10) 



with 



a 



p.q.r.s = rs - pq\ bp_,i,-,.s = 6pqrs{r + s-p-q)+4a [rs{r + s) ~ pq{p + q)] +9a{rs - pq); 
'p.7/.T.i = rs-pq; bp,gj,-s = 6pqTs{T + s -p -q) +4a[Ts{T + s) -pq{p + q)] +9a{Ts-pq); 



and ap^q,r,sb-p;qj-s ^ a-p;cj;frsbp^q.r/, p, q, r, s; p, q, r, s ai-e integers. Then 
Proof. According to the proof of Theorem 12. 2 1 we have 



Hcof 



J 6pqrs{r + s-p-q)+4a[rs{r + s)-pq{p + q)]+9a{rs-pq) 3,^/4 

36{rs — pq) 



z^ + ff{z*) u{0}) 



and 



=^(3Mwj«)(®) 



P,1,' 

= i-icor 



6pqrs{r + s-p-q)+4a[rs{r + s)-pq{p + q)]+9a{rs-pq) 3 , ^. 4-, ^r ^ 

1 + L/— X Z^ + ff{z ) u{co). 

io(rs — pq) 



Then there exists 

and by the similar way to the proof of Lemma |273] we get 
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n 

For the right Riemann-Liouville fractional derivative, denote that 

«A««W--lf;^f«(x+(fe-/,)/z), (2.11) 

" A:=0 

where q^ is defined by ( I2.2l i and p an integer Using the same way as Theorems 12. 1112.31 we can obtain 
the following results. In particular, the coefficients in (12.121 ) are completely the same as the ones in 
( 12.61 ): the coefficients in ( 12.131 ) the same as the ones in ( 12.71 ); and the coefficients in ( 12.141 ) the same as 
the ones in (12.81). 



Theorem 2.4 (Second order approximations for right Riemann-Liouville derivative) Let u, xD"'^^u{x) 
with a e (1,2) and their Fourier transforms belong to Li (M), and denote that 

2RAp^gU{x) = WpRApu{x)+WqRA"u{x), (2.12) 

then 

Theorem 2.5 (Third order approximations for right Riemann-Liouville derivative) Let u, ^D"'^^u{x) 
with a e (1,2) and their Fourier transforms belong to Li (M), and denote that 

3RAp.y.r..Ax) == Wp^g2RA",^u{x)+Wr,s2RA%u{x), (2.13) 

then 

,D2m(x) = 3«A,%,,m(x) + ^(/i3). 

Theorem 2.6 (Fourth order approximations for right Riemann-Liouville derivative) Let u, xDZ'^'^u{x) 
with a E (1,2) and their Fourier transforms belong to Li (M), and denote that 

4ffA"^ ,,,P,^,5=,jM(x) = Wp,^,r.,3M"^,,-,M(jc) + Wp,5,7,I3/fAf,5,7jM(-«), (2-14) 

then 

All the above schemes are applicable to bounded domain, say, (x^, xr), after performing zero ex- 
tensions to the functions considered. Let u{x) be the zero extended function from the bounded domain 
(jci, Xr), and satisfy the requirements of the above corresponding theorems ( Theorems 12. 1112. 6V Denot- 
ing 

1 1^1+'' 
LA^pUix)^- ^ q^u{x^ik-p)h), (2.15) 

" k=0 

then 

,,Z)Xx) = iA^x) + ff{h), p^O, 

~ r, (2.16) 

;,,DXx) = lA,Xx) + ^(/z2), p = 0- 
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xi^Oyix) = 2lA"^m(x) + ff{h^), where 2lA"^m(x) = WpLA'j^u{x)+WgLA"u{xy, (2.17) 

xlD^u{x) = iLAp^^.r.A^) + ^{h^)^ where 3iA;]^^,,M(x) = Wp,^2lA"^m(x) + W;-,,2lA"^m(x); (2.18) 

and 

,,D>(x) = 4iA«^ ,.,p_^,,,«(x) + ^(/i^), (2.19) 

where 4iA^«^,,j,^,,jM(x) = Wp,^,,.,3lA,"^^,.,m(x) + Wp,5,r,j3LAf5,7,s«(-^)- 

Denoting xi =xl + ih, i = —m, . . . ,0,1,. . . ,Nx — i,Nx, ...,Nx + m, and h — {xR — xl)/Nx being the 
uniform space stepsize, it can be noted that 

u{xi)—Q, for / = — m,— m+ 1,...,0 and i =Nx,Nx+l,---,Nx + in, 

where 

m = max{abs{p,q,r,s,p,q,r,s)). (2.20) 

Then the approximation operator of (12.15) can be described as 

^ 1 i+p 1 i+m 1 i+m 

^""(•^') = I^ E ^""(■^'-*:+p) = I^ E ?"+p-mM(-^<-A:+m) = T^ E ^"+p-m«(-^<-A;+m), (2.21) 
i:=0 k=m—p k=Q 



where qf,„_,„ = 0, when k + p — m <Q, and /9 is an integer Then 

1 



HP 



,,DX^,) = lA«m(^,) + ^(/i) = T^ L ??+p-™«(-^'-*+m) + ^(h), p^O, 



k=0 

i+m 



-^ 1 I+m 

x.Oyixi) = LA"pU{xi) + ^(/i2) = — J2 ?f+;,_^M(^,-^+m) + ^{h') , /^ = 0; 



(2.22) 



uc 



~ 1 '+™ 

,Dy{xi) = 2lA«^«(x,) + ^(/z2) ^ — Y, {yvpqt+p-,n + vv,?:^+,_,„)«(x,_,+„) + e(h\ (2.23) 



" yt=0 






1 i+m 
*:=0 



+ ^(/i2); 

1 i-\-m 

x^D^uixi) = ALA%^r,s.p,W.A''i) + ^('^^) = I^ E ^"(^i-t+m) + ^(^^). (2-25) 

" /t=0 



where 



+ Wp,^,,vvWr..s w.v?"+j_„, + Wp.f .7,jWp,;jWp^"+p_„, + Wp,5,r,jWp,5W^^f+^_^ (2.26) 

+ WpJi^-sWj-sWTq'j^^j_^ + W-p^j;sWrrsWjq'j^+-,_„,. 

Taking f/ = [m(xi),m(x2), • • • ,m(-Xa(j-i)]^, then (12.21b can be rewritten as the matrix form 



^ -h 



4 



-k 



lOofEU 

where 



A" 



<+2 



q"-2 



^p- 1 ^0 

^p '/p-l 



Hp+n-3 

a a 

yp+n-l qp+n-i 



a 



?o" 






?o" 



a a a 

■ 1p+l <ip 1p~l 

€-2 ■■■ <+2 <+l <. 



and /? is an integer and q" — 0, when A: < 0. From (I2.23b -( l2.25l l we obtain 

iLA-pg^JJ = T^^p,q,r,s^: ^p.q.r.s — ^P-.g^P.'l + Wr,sAr.s', 

Xa r r <a Tr ^a ^a -i^u, A'^ 

4L-f^p^qj\s.p,q,r.s'^ — .aP-^-T^l^^JJi^ ^p.,g.,r,s,p.;g.r.]i ^ '^P.'i.r.s-^p^q^r.s ^ '^P.g.r.s-^p,?/.?,!- 

Similarly, for the right Riemann-Liouville derivative, taking 



(2.28) 



(2.29) 

(2.30) 
(2.31) 



rA^m(x)==— ^ q^u{x+{k~p)h), 



k=Q 



then there exists 






h'- 



k=0 



k=m—p 



" k=Q 



where q'^^„_„, — 0, when k + p — m<Q, and p is an integer And the fourth order approximation is 



xD^A^i) = 4RA'j',^q,r,s,M,7A''') + ^C""^) 



1 Ni — i+m 



Y, 9ku{xi+k- 



k=0 



where cp" is defined by (12.261 1. and the matrices forms are 



(2.32) 



^KU = T^B-U, B- = (A 



a\T. 



2RKaU = TwK.aU, B" = W.B" + W,B"; 



P:1 



\a P-g 



p.q - '"P^p ^ '"q^q ■■ 



3RAp^g,r.,U - J^B"^,.^U, -Bp.^.r^,. - Wp.qBp^q + Wr.sBr/, 
'iR^-p q^^-pljfjU ~ TtB p^^ ^-piifsU , Bp i^^.^-p-qj:j = W p^q^r,sB p q y s + W-pqjja-p.qy:-g. 



(2.33) 



-v 



-h 



-4- 



-k 



llof|24] 



Remark 2. 1 When /? = 0, then A" in ( I2.28l l reduces to the lower triangular matrix, and it can be easily 
checked that all the eigenvalues of A" are greater than one; in fact, from Lemma 12.11 it can be noted 

that A(AS:) = (I) , with a G (1,2). This is the reason that the scheme for time dependent problem 
is unstable when directly using the second order Lubich formula with a G (1,2) to discretize space 
fractional derivative. 



2.2 Effective fourth order discretization for space fractional derivatives 

This subsection focuses on how to choose the parameters p,q,r,s,]5,q,T,s such that all the eigenvalues 
of the matrix A" (orA" ,.,. orA" .^----) have negative real parts; this means that the corresponding 
schemes work for space fractional derivatives. Since B" B" j.^, and B" rspqri i^' respectively, the 
transpose of A" , A" ^^, and A" .^. ^^-j, we don't need to discuss them separately. 

Definition 2.7 dOuarteroni et a/.ll2007i p. 27) A matrix A G M"^" is positive definite in M" if (Ajt:,jt:) > 
0, Vjc g M", X 7^ 0. 

Lemma 2.4 JOuarteroni et a/.L 120071 p. 28) A real matrix A of order n is positive definite if and only if 



its symmetric part H - 



, A+A' 



is positive definite. Let H G M"^" be symmetric, then H is positive definite 



if and only if the eigenvalues of H are positive. 

Lemma 2.5 JOuarteronief a/.[l2007l p. 184) If A G C"''", let H = ^^^ be the hermitian part of A, then 
for any eigenvalue A of A, the real part 9{(A(A)) satisfies 

K,„{H) ^^{X[A)) ^Ka,{H), 

where Xmin{H) and X,„ax{H) are the minimum and maximum of the eigenvalues of H, respectively. 



Definition 2.8 dChan & Jinl l2007l p. 13) Let n x n ToepHtz matrix T„ be of the following form: 



to 



t-i 
to 



t-i 

tl to 

'■■ tl 

tn-2 

tn—\ ^«— 2 



i.e., tij — ti^j and T„ is constant along its diagonals. Assume that the diagonals {f/t}^=i„ i i are the 
Fourier coefficients of a function /, i.e.. 



t2-n 


t\-n 




t2-n 




t-\ 


t\ 


fo . 



1 



tk = ^J_Jix)e-"''dx, 

then the function / is called the generating function of T„ . 

Lemma 2.6 (IChan & JinL 120071 p. 13-15) (Grenander-Szego theorem) Let r„ be given by above matrix 
with a generating function /, where / is a 2;r-periodic continuous real- valued functions defined on 
[— TT, n]. Let X,„m{Tn) and X,„ax{T„) denote the smallest and largest eigenvalues of T„, respectively. Then 
we have 

J win ^ ^miny^n) ^ ^maxy^n) ^ fmaxj 



-i- 



-h 



4 



-k 
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where /„„■„ and f,„ax is the minimum and maximum values of f{x), respectively. Moreover, if /„,„ < 
fnmx, then all eigenvalues of T„ satisfies 

Jmin ^ ^\^n) "^ Jmaxj 

for all n > 0; In particular, if /„„„ > 0, then T„ is positive definite. 



Theorem 2.9 (Effective second order schemes) Let A"„ be given in ( 12.29b and 1 < a < 2. Then any 
eigenvalue A of A" satisfies 

9^(A(A«J)<0 for {p,q)^(l,q), \q\^2, 
moreover, the matrices A" and (A" )^ are negative definite. 
Proof. 
(1) For ip,q) ^ {l,q), q ^ -2, we have A",, = ^{qAf -A«), and 



with 





r <p^ 




C 






.a _ 














C-2 


K-2 




0« 


K 
K. 


tf = 


r 11? 
) 1-1 


_.,a 

1k+q-i 

y-1 


k> -q. 


-q, 



The generating functions of A" and (A" )^ are 



fA?,ix) = L tf ^'^'^''^ ^nd /(^„^).(x) = ^ tf . 



-i(/t-l).t 



/t=0 



<:=0 






is the generating function of 



respectively. Taking i/p,^ = :l££±il££L, then/p,^(a,x) : 
i/p,^. Since fji^a^ (x) and /mh y [x] are mutually conjugated, then fpg{a,x) is a 27r -periodic continuous 
real-valued functions defined on [— ;r,;r]. Moreover, fp^g{a,x) is an even function, so we just need to 
consider its principal value on [0, n]. Next, we prove fp^q{a,x) ^ 0. Rephrasing the generating function 
leads to 

^ \k=Vi k=0 / 

I ^1 \ k=Q k=Q k=0 k=0 



a „—ikx 



2{q- 

1 

2(<7-l) 
1 



^e-«(l _e")« 1 + -(1 -e«) +qe'^{\ 



-,x^ai 1 + 1(1 



2(^-1) 



-"'■^(1 -e'-^)« ( l + -{l-e'-')] +e'i''{l-e-'-Y M +1(1 -g-'^) 



-h 



-h 



-4- 



-k 



Because of 



where 



9 = larctan- 
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(1 -e±-)« = (2,««^)%±'«(i-?), ^+1(1 -e±'- 



2s;n 7 



cos| + wl +3sin^^ 



= {l+3sm'-ye^'"^i-'' 



e[0,7r/2], 



then, for ^ ^ —2, there exists 

1 / x\ « / 
fp^y{a,x) = —— \2sin-j \l+3sm 



2. 



^COi 



(a(x-| 



0) —X 1 — cos 



(2) For (p,^) = (1,^), <7 > 2, we have A«^ = -^(<?Af -A«) and 



( a(x— — — 9) — qx 






€ 



'9+2 









C 



Ci 



c 



c 



€-2 






<+l -^.^ Cl 



C 



K-2 ••• ^% Ci <J 



with 



tf = 



-^i^, O^k^q-2, 



-^3j , k>q-2. 



The generating functions of A"g and (A"^)^ are 

00 00 

/^%W = L tf ^'■''-'''■'^ and /(A«,)K^) = E 0fe-'('-^)^ 
respectively. Denoting 



yt=0 



i=0 



^P,'? - 



^p,<:/ + \^p..q) 



(2.34) 



then fp,q{a,x) = —^ y'^^ is th^ generating function of Hpq. By the similar way, for q^2. 



there exists 






Q'coi f a(x 0) — x\ —COS I a{x 0) 



- qx 



-^■ 



-h 



4 



-k 
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Fig. 1. /(a,x) for (p,^) = (1,-2) 



Fig. 2. f{a,x) ^ for {p,q) = (1.2) 



It can be noted that fp,g{a,x) has the same form when q ^ —2 and q ^ 2, p — 1. And we can check 
that, for {p,q) = (!,<?), |^| ^ 2, there exists (see Figs.[T]|2]l 



fp,q{a,x) = 



qcos I a{x Q) — x) —cos I a(x Q) — qx 



(2.35) 



s^O. 



Since fp^g{a,x) is not identically zero for any given a £ (1,2), from Lemma l276l it implies that 
X{Hp,g) < and Hp,^ is negative definite. Then we get 9{(A(A" )) < from Lemma |231 and the 
matrices A"^^ and (A"^)^ are negative definite by Lemma l24l D 



Theorem 2. 10 (Effective third order schemes) Let A"^ ,. ^ with 1 < a < 2 be given in ( I2.30l i. Then any 
eigenvalue X of A" ,. ^ satisfies 

9^(^(^",,J)<0 for {p,q,r,s)^{l,q,l,s), |^| > 2, |5| > 2, and ^. < 0; 

moreover, the matrices A" ,. ^ and (A" .)^ are negative definite. 
Proof. Taking 



4" 

^p.q.r.s — 



V^p.q.i 



— Wp.qHp^q + WrjHr^s 



where Hpq and i/r,s are defined by ( |2.34t . then 



fp.q,rA<^,x) = Wp^q fp^g{a,x) +Wr,sfr.s{CC,x) 



(2.36) 



(2.37) 



is the generating function of i/p,^,,-,.,, where fp,g{a,x) and /„■(«, x) are given by (12.35b . Since |^| ^ 2, 



\s\ ^ 2, and ^.s < 0, we can check that Wp^g — wi^q — y _ s 
(l235T l and (lOTI i. we get fp^q^,.A<^,x) «: 0. 



> 0, Wr 



Again, from Lemmas 12.4112. 61 the desired results are obtained. 



n 



-h 



-h 



A- -k 
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Theorem 2.1 1 (Effective fourth order schemes) Let A" rspqri ^i* 1 < a < 2 be given in ( 12.311 ), 
where {p,q,r,s,p,q,T,s) = (1,2,1,-2,1,^,1,1), \q\ ^ 2, \s\ > 2, (q,s) ^ (2,-2) and qs < 0. Then any 
eigenvalue X of A"^ ^^.p^- j satisfies 

and the matrices ^"^.^.^,^.5.7.1 and (A"^^ ,. ^ -p^- j)^ are negative definite. 
Moreover, if {p,q,r,s,p,q,T,s) takes the following values 

(/?,^,r,s,7?,^,r,S) = (1,2,1,0,1,2,1,-2), 
{p,q,r,s,p,q,T,-s) = (1,2,1,0,1,-1,1,-2), 
{p,q,r,s,p,q,T,-s) = (1,2,1,-1,1,2,1,-2), 
{p,q,r,s,p,q,r,s) = (1,2,1,-1,1,-1,1,-2), 
{p,q,r,s,p,q,r,s) = (1,0,1,-1,1,2,1,-2), 
{P,q,r,s,p,q,r,'s) = (1,0,1,-2,1,2,1,-2), 
{p,q,r,s,p,q,r,s) = (1,-1,1,-2,1,2,1,-2), 

then 5R(A(A«^ ^,^^--)) < and the matrices A«^ ,.^.---- and (A«^ ^,,----)^ are negative definite. 
Proof. By the similar way to the proofs of Theorems 12. 9 1 and 12.101 we obtain the desired results. D 

3. Application to the space fractional diffusion equations: the one dimensional case of ( 11.2) and 
dO) itself 

We use two subsections to derive the full discretization of (11. 2t . First, we present the scheme for the 
one dimensional case of (11.2b . The second subsection detailedly provides the full discrete scheme of the 
two-dimensional fractional diffusion equation (ll.2l l with variable coefficients. 

3.1 Numerical scheme for ID 

In this subsection, we consider the one-dimensional case of (|1.2| i with variable coefficients, namely, 

du{x,t) 



dt 



^ d+{xl„^Dy{x,t) + d-ixy^D"u{x,t) + f{x,t). (3.1) 



In the time direction, we use the Crank-Nicolson scheme. The fourth order left fractional approximation 
operator (I2.25t , and right fractional approximation operator (I2.32l i are respectively used to discretize the 
left Riemann-Liouville fractional derivative, and right Riemann-Liouville fractional derivative. 

Let the mesh points x, = xi + ih, i = —m, ...,0,l,...,Nx — l,Nx,---,Nx + m, where m is defined 
by (12.20b and f„ = nT, ^ « ^ A^(, where h = {xr — xi)/Nx, i = T /Nt, i.e., h is the uniform space 
stepsize and T the time steplength. Taking m" as the approximated value of u{xi,t„) and £/+,,■ = d+{xi), 

d^i — d-{xi), f" = f{xi,tn+\/2), where f„+i « = (f,, +r„+i)/2. Then, Eq. (13.1b can be rewritten as 



u{xi,tn+l) - u{xi,t„) 1 



d+.i 4LApg,-sjj^gj^jU{xi , tn+ \)+d+^i 4LApyr,s.j},7j,r.l'^{^' ; ^n ) 

+ d-jARA-pgr.s.p.Ji.r.l'^i^i^fn+l) + d-^ix^ ARApgj-^s^p.jj.T.l'^iXiJn) 
+ f{xi,t„+,l2) + ff{T^ + h''). 



(3.2) 



^ -h 



4 
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Multiplying ( I3.2l i by T, we have the following equation 



with 



1 + X ( "+,i4L-Ap^ r^,^^^^jj + a_.i4KAp_^ ,.5 p_^j=j I u[Xi,t„) + T/(Xi,f„_|_i/2) +-'^' 






|7;''+l|<?TfT2 + /z4) 



Therefore, the full discretization of ( 13.3) has the following form 

J \^+,i'iL^p,q,r,s,'p.;q.T.s ~^ '^-,i4-R^p,q,r.s.;p.;q.T.s j ^i 



and it can be rewritten as 



„n+l 



J i+;n J Ny — i+m 

^^ Li'Pk "i-k+m +T^ L ^k "i+k-ii 



k=0 



k=Q 



J i+m J Ny-i+m 



h 



h" 






«+l/2 



4=0 " A-=0 

For the convenience of implementation, we use the matrix form of the grid functions 

J J" — ],,",," „n iT 77«+l/2 _ [^"+1/2 ^n+1/2 Ji+l/2nT 

U — [Mi,M2,---,«Af,-lJ , r ' —[J\ ,J2 '■■■'/Afj-l J ' 

therefore, the finite difference scheme (13.6) can be recast as 



I~^{D+Aa+D^A',} 



U 



n+l 



I+ — iD+Aa+D^Ai)\u'' + TF"+'/\ 



where Ace = A" .^.^^j=j is defined by (12.31b . and 



D+ = 



U.2 



' + .N,-l 



, D = 



-,N,-l, 



(3.3) 



(3.4) 



(3.5) 



(3.6) 



(3.7) 



(3.8) 



3.2 Numerical scheme for 2D 

We now examine the full discretization scheme of (11.2b . For effectively performing the theoretical 
analysis, we suppose d+{x) — d+{x,y), d-{x) — d-{x,y), and e+{y) = e+{x,y), e^{y) — e^{x,y). 

Analogously we still use the Crank-Nicolson scheme to do the discretization in time direction. 
Let the mesh points x,- — xl + ih, i = —m, ...,Q,l,...,Nx — l,Nx,---,Nx + in, and yj = yL + j^y, j = 
—m, ...,Q,l,...,Ny — l,Ny,...,Ny+m, where m is given in ( 12.20b . t„ — m, Q^n^Nt, and Ax = [xr — 
xl)/Nx, Ay = {yR-yL)/Ny, t = T/N,; and <i+.,- = d+{xi,yj), d-j = d-{xi,yj), and e+j = e+{xi,yj). 



-v 



-h 



-4- 
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e_ j = e_(x,-,y,). Taking u" ■ as the approximated value of u{xi,yj,t„) and /". ' = f{xi,yj,tn+i/2}' 
where f„+i/2 = {t,, +t„+i)/2. Then, Eq. ( 11.2b can be rewritten as 



+ Tf{x„yj,t„+i/2)+Rij\ 



with 



Then, the resulting discretization of (|3.9l l has the following form 



(3.9) 
(3.10) 



1 ~ 2 ('^+'''*^P,'?.'-.*:M.r.i "•" '^->'4K^p,<?,r,i,p,9,7J + ^+:./4i^p,(?,r,i.75,5.7,J + ^- J '^R^p.q.r.sJjJjJ-rs 



/'+^ 

^',; 






We further define 

thus Eq. ( 13.111 ) can be rewritten as 

The perturbation equation of (13.12) is of the form 

(l - |5„,.) (l - l5p,) <+' = (l + 15„..) (l + l5p,) ««,. + T/;:^ 
Comparing (13.13b with ( 13.12b . the splitting term is given by 



,s.p,q,r.s 



'J 



(3.11) 



(3.12) 



(3.13) 



-5„,.5^,,«r-<j) 



since {u"j — u" ■) is an ^{t) term, it implies that this perturbation contributes an ^(t ) error compo- 
nent. 

The system of equations defined by ( 13.131 ) can be solved by the following schemes. 

PR-ADI scheme IPeaceman & Rachford ( 1955) 1: 



(l-|5„,)<, = (l + |^,)<, + |/^'/^; 



(3.14) 



-i- 
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D-ADI scheme llDoudJ ( Il955h l: 

(i - 15„,,) uij = (i + 15„, + T5p,,) <,. + T/;'7^/^ 



l-o^/^'vXT =<j-^hAr 



(3.15) 

(3.16) 
(3.17) 



Take 



IT" — III" u" ii" u" /y" ii" u" ii" ij" 1^ 

TT'n r fn ftl rn fn fn rn fn rn fu l T 

•^ — [/l.li/2,li--- '/«,--!, 17/L27/2,2i---'/A'^-1,2j---i/1,Wv-1'/2.A'v-1 J- ■■i/Wv-LA'y-lJ ) 

and denote 



<. 



2{Ax) 

T 

2(4^ 



-[(/®D+)(/®A„) + (/®D_)(/®A^)]=— -— /®(D+A„+D_A^ 



(£+ /) (A/3 (g>I) + (£_ (g) /) (A^ ® /) 



2(zix) 

T 
2(4^ 



£+A«+£_Ai ®/ 



(3.18) 



where / denotes the unit matrix and the symbol (g) the Kronecker product [see, iLaua (120051) 1. and Aa = 
A" ■p-^--^,Ap —Ar^-p-j-^ are defined by ( 12.311 1. The matrices D+ andD_ are defined by ( I3.8l l, and 



'+,1 



e+,2 



£■_ 



Therefore, the finite difference scheme (13.13) has the following form 



(3.19) 



(3.20) 



Remark 3.1 T he schemes (I3.14li-(l3.15l l and ( I3.16I )-( I3.17I ) are equivalent, since both of them come 
from (ITTJI i. see JDeng & Chenl(l2013l) l. 



4. Convergence and Stability Analysis 

In this section, we theoretically prove that the difference scheme is unconditionally stable and 4th order 
convergent in space directions and 2nd order convergent in time direction. In the following, the matrices 
D+, D and E+, £_ are defined by (13.81 1 and (|3.19t , respectively. 

Lemma 4.1 (lLaubll2005[ p. 140) Let A e M'"^", B e R''^', C e R"^'\ and D e R""". Then 



{A(g)B){C(g)D)^AC(E)BD {eR""'''"). 
Moreover, for all A and B, (A B)^ = A^ ® B^. 
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Lemma 4.2 ( lLaubll2005[ p. 141) Let A e W" have eigenvalues {l,}^^i and B e R'"""" have eigenval- 



ues {/XjI'J'^j. Then the mn eigenvalues of A (g)B are 

Ai^l , . . . , Ai;U„,, ^2^1 , . . . , A2/X„,, . . .,X„Hi . . .,XnjJ.,„- 

Theorem 4.1 Let the matrix A „ = A" .^-^-.^ be defined by (12.31b andD_ = KaD+, where Ka is any 
given nonnegative constant. Then we have 5R (A {D+{Aa + K"aA^))) < 0. 

froo/ Since 

D7 [D+(Aa + (C„A^)] Dl ^ DliAa + KaAl)Dl, 

1 1 

it means that D+{Aa + K'a^a) ^'i'^ 0+(^a + '^aA^)D\ are similar From Theorem 12. Ill we know Aa 
and A„ are negative definite, and thanks to Definition l2.7l it implies that 



Dl{Aa + KaAl)Dlx,x\ = i{Aa + KaAl)Dlx,Dlx\ <0, VxGM",x^O, 

i.e., the matrix A := D\{Aa + KaA^)D\ is negative definite. From Lemma l24l H — ^^^ is negative 
definite and ^ax{H) < 0; and according to Lemma l231 we obtain 5R(A(A)) ^ Xmax{H) < 0. Therefore, 

9^(A(D+(Aa + fc„A^)))=5R(A(A))<0. D 

Theorem 4.2 Let ^ and £/y be defined by ( 13.181 ) and D_ = KaD+, £_ = «•/}£'+, where Ka and JC^ are 
any given nonnegative constants. Then we have 3i (A (is^)) < and 3i (A {^/y}) < 0. 

Proof. From (13.18b . there exists 

By Theoremim weget5R(A(D+(Aa + K-„Aj))) <0 and3i(x (E+iA^ + KiiAJ)]) < 0. Then, ac- 
cording to Lemma|43] it implies that 9^ (A(j/x)) < and 5R {X{£/y)) < 0. D 

Remark 4.1 If taking fCa = kj3 = 0, then Eq. (11.21 ) becomes the one-sided fractional diffusion equation; 
and Ka — Kn — l, Eq. ( 11.2b reduces to the space-Riesz fractional diffusion equation. 

4. 1 Stability and Convergence for ID 

Theorem 4.3 Let D = KaD+, then the difference scheme (13.7b with a G (1,2) is unconditionally 
stable. 

Proof. Let S]' (/ = 1 , 2, . . . , A^^ — 1 ; n = 0, 1 , . . . , M) be the approximate solution of u", which is the exact 
solution of the difference scheme (13.7b . Putting ef = S^' — m", and denoting e" = [£",£2,... ,ej^ _J, then 
from (13.7b we obtain the following perturbation equation 

(/-A)£"+i = (/+A)£", 

i.e., 

£"+i = (/-A)-i(/ + A)£", 



^ -h 
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with 

A = ^D+{Aa + KaAl). (4.1) 

Denoting A as an eigenvalue of the matrix A, then from Theorem 14. II we get 3i{X{A)) < 0. Note that 
A is an eigenvalue of the matrix A if and only if 1 — A is an eigenvalue of the matrix I — A, if and only 
if (1 -A)"' (1+ A) is an eigenvalue of the matrix {I -A)~^{I +A). Since 9^(A(A)) < 0, it implies that 
1(1— A)^^(l+A)| < 1. Thus, the spectral radius of the matrix (/— A)^'(/+A) is less than 1, hence the 
scheme (13.7b is unconditionally stable. 

D 

Theorem 4.4 Let u(xi,t„) be the exact solution of ( 13.11 ) with a G (1,2), m" the solution of the finite 
difference scheme (13.7b . and D = KaD^, then there is a positive constant C such that 

||M(x,-,f„)-<||2S^C(T2 + /;4), /=l,2,...,A?,-l;n = 0,l,...,M. 

Proof. Denoting e" — u{xi,t„) —uf, ande" = [e'jjgj,- ■ ■ ,e5(, j]^. Subtracting (13.2b from (13.7b and using 
e" = 0, we obtain 

{I-A)e"+^ = {I + A)e"+R"+\ 

where A is defined by ( 14.11 ). and R" ~ [/?" ,/?2, . . . jR"^ -iV ■ The above equation can be rewritten as 

e"+i = (/-A)-i(/+A)e" + (/-A)-i7?"+i. 

Similar to the proof of Theorem 4.2 of toeng & ChenI (|2013|) 1. we have that ||(/-A)-'(/+A)||2 and 

||(/-A)-i||2 are less than 1. Then, using |7?^+'| < ?'t(t2 +/j4) in llHl, we obtain 

||."||2<||(/-A)-'(/ + A)||2.||."-'||2 + ||(/-A)-'||2.|r| 

^\\e"-'\\2 + \R"\^"Y,\R'+'\^cix^ + h% 

k=0 

n 

4.2 Stability and Convergence for 2D 

Theorem 4.5 Let D = KaD+ and £_ = KpE+, then the difference scheme ( 13.20b with 1 < a, j3 < 2 
is unconditionally stable. 

Proof. Let «)' • {i = 1,2, .. . ,Nx — l',j = 1,2, .. . ,A',^ ~ l;n = 0, 1, . . . ,Nt) be the approximate solution 
of uf j, which is the exact solution of the difference scheme (13.20b . Taking e" =ufj~ m" , then from 
( 13. 20b we obtain the following perturbation equation 

[I ~ £^,)ii ~ s^y)e"+^ = {I + s^,){i + s^,)e\ (4.2) 

where s^x and jz/y are given by (13.18b . and 

en r„n „« „« „n „n „n „h „n „« iT 

— [ti,i,fc24,---,fcArj_i,i,ti,2it2^2J---'tAr^-l,2'---'ti^;V,— ljt2,w,,-l,---,fcAfj-l,W,,-lJ ■ 

Then Eq. ( 14.21 ) can be rewritten as 

£"+' = [I - £/,)-\i ~ s^xY'ii + s^x){i + <)e" ■ (4.3) 
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According to Lemma l4n and ( 13.181 1. it is easy to check that ^/^ and £/y commute, i.e.. 



s^.j^y = £/,.£/, = ^ ^^^ p {s+Ap + E^Aj'^ ® {D+Aa + D^Al) . (4.4) 

Then Eq. ( 14.31 1 has the following form 

£"+' = {I - £/,)-\l + £^:,){I - £^^-\l + J^-)e" . (4.5) 

Form Theorem l4.2l we have 5K (A(£^)) < and 3i (A(^.)) < 0. Similar to the proof of the Theorem 
14.31 the spectral radius of the matrix (/ — ^) ^ ' (/ + £^x) and (/ — J2^.) ^ ' (/ + £^.) are less than 1 . Then 
the difference scheme (13.20b is unconditionally stable. 

D 

Theorem 4.6 Let u{xi,yj,t„) be the exact solution of ( |1.2t with 1 < a,l5 < 2, u" ■ the solution of the 
finite difference scheme (13.20b . and D = KaD+ and £_ = KpE+, then there is a positive constant C 
such that 

\\u{xi,yj,t„) - ulj\\2 < C(t2 + (Zix)4 + (Ziy)^), 

with / = 1,2, . . . ,A^, - l;y = 1,2, . . . ,Af>. - 1; n = 0, 1, . . . ,A^,. 

Proof. Taking e" ■ = u{xi,yj,tn) — m" , and subtracting (13.9b from (13.20b . we obtain 

{I ~ .s/,){I - .s/y)t"+' = {I + £/,){! + £/y)e" + R"+\ (4.6) 

where jz^ and £/y are given in ( 13.18b . and 



1, 1' ^2,1 :■■■ 1^^,-1,1 ;^l,2:^2,2:---i^W,-l, 2'- ■■ 1^1, A'y-l'^2,A',,-b---'^A',— l,A',,-lJ J 

P" P" P" 

1,2; • • • j"l,Wy-l;"2,Wj,-l J ■ ■ ■ :"A(r-l,A^y-lJ ; 



Rn I nil pii TfU pi p" p« pn p" p« 1^ 

— ["l,l;"2,lJ---:"Arr-l,l'"l,2j"2.2i---'"A',-l,2!--- J"l,^'y-l^"2,M,-lJ■■■:"^ 



and \R'l+\ < ?t(t2 + (zix)'* + [Ayf) is given in dlTol ). 

From (14.4b . ^ and ^, commute, then Eq. ( 14.61 ) can be rewritten as 



(/-^J-i(/ + ^,)(/~^y)-i(/ + <)e" + (/-^,)-i(/-^y)-^R" 



Again, similar to the proof of Theorem 4.2 of lIDeng & ChenI (1201 3h l. we know that ||(/ — a^v) ' [I - 

^v)||2and 11(7 — jz/y)^ lb are less than 1, where V=x,y. Then there exists 

k=0 



D 



5. Numerical results 



In this section, we numerically verify the above theoretical results including convergence rates and 
numerical stability. And the loo norm is used to measure the numerical errors. 
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5.1 Numerical results for ID 

Consider the one-dimensional fractional diffusion equation ( 13.1b in the domain 0<x<2, 0<f^ 1, 
with the variable coefficients d+ [x) —x"',d^ [x) — 2x", and the forcing function 

r(9) 



f{x,t)^cos(t + l)x^{2-xf-x'^sin{t + l) 



S-HN 



r(9-a) 



(;c«-«+2(2-jc)''-«) 



res -a 



.(x 



1-a 



-2(2 -x)'-") +24 



r{i) 

r{l-a) 



6-a\ 



(x''-"+2(2-x)*'-«) 



■32- 



r(6) 



-s^- 



,5-a 



-2(2 -jc)^-«) + 16^^^^(x4-« +2(2 -x)4-«) 



r{e-a) ' ' ' r(5-a) 

and the initial condition u{x,Q) ~ j'/«(l)x'*(2 — x)^, the boundary conditions M(0,f) == M(l,f) = 0, and 
the exact solution of the equation is u{x,t) — sin{t + \)x'^{2 —x^. 



Table 1. The maximum en'ors and convergent orders for the scheme 13. 7t of the one-dimensional fractional diffusion equation 
(3lT)att=l andT = /i2. 



{p,q,r,s,p,q,r,s) 


h 


a = 1.1 


Rate 


a = 1.9 


Rate 




1/10 


4.7842e-03 




5.8264e-03 




(1,2,1,0,1,2,1,-2) 


1/20 


2.5436e-04 


4.2333 


5.9999e-04 


3.2796 




1/40 


1.9662e-05 


3.6934 


4.6242e-05 


3.6977 




1/60 


4.1748e-06 


3.8218 


9.7725e-06 


3.8334 




1/10 


8.5475e-03 




5.5003e-03 




(1,2,1,-3,1,2,1,-2) 


1/20 


4.9722e-04 


4.1035 


5.7476e-04 


3.2585 




1/40 


3.9559e-05 


3.6518 


4.4490e-05 


3.6914 




1/60 


8.6604e-06 


3.7464 


9.4148e-06 


3.8301 



Table \T\ shows the maximum errors, at time t — 1 with T = h^, the numerical results confirm the 
convergence with the global truncation error €?(t^ + /z^). 

5.2 Numerical results for 2D 

Consider the two-dimensional fractional diffusion equation ( 11.2b . where 0<x<2, 0<3'<2, and 
<t ^ 1, with the variable coefficients d+{x,y) —x", d-{x,y) —2x", ande+(x,y) =}'", e^{x,y) ~2yP, 
and the initial condition m(ji:,3;,0) = sin{l)x'^{2 — x)y (2 — y)"* with the zero boundary conditions, and 
the exact solution of the equation is 

u{x,y,t) = sin{t + 1)^(2 -x)V(2 -y)4. 

From the above conditions, it is easy to get the forcing function f{x,y,t). 

Table |2] displays the maximum errors of the scheme (13.20) . and confirms the desired convergence 
with the global truncation error ^(t^ + (Ax)^ + {Ay)^). 

6. Conclusions 

Based on the Lubich's operators, this work provides a new idea to obtain the high order discretization 
schemes for space fractional derivative. We obtain the effective difference operators with 2nd order. 
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Table 2. The maximum en'ors and convergent orders for the scheme 13. 201 of the two-dimensional fractional diffusion equation 
fTH at t=l and T = {Axf = {Ayf. 



{p,q,r,s,p,q,r,s) 


Ax 


a = /3 = l.l 


Rate 


a= 1.8,jS = 1.9 


Rate 




1/10 


8.6154e-03 




6.521 le-03 




(1,2,1,0,1,2,1,-2) 


1/20 


5.4115e-04 


3.9928 


4.4802e-04 


3.8635 




1/30 


1.2626e-04 


3.5894 


8.8416e-05 


4.0023 




1/40 


4.3328e-05 


3.7177 


2.7791e-05 


4.0229 




1/10 


l.OllOe-02 




6.6368e-03 




(1,2,1,-3,1,2,1,-2) 


1/20 


6.388 le-04 


3.9842 


4.547 le-04 


3.8675 




1/30 


1.4363e-04 


3.6806 


8.9704e-05 


4.0032 




1/40 


4.8431e-05 


3.7788 


2.8199e-05 


4.0226 



3rd order, and 4th order accuracy, called WSLD operators. For further checking the efficiency of the 
high order schemes, we apply the 4th order scheme to solve the space fractional diffusion equation with 
variable coefficients; and the detailed theoretical analysis and numerical verifications are presented. 
Hopefully, the higher order (5th order, 6th order, etc.) schemes can be obtained by following the idea 
given in this paper In fact, for any fixed convergent order, the obtained difference operators are a class 
of difference operators, not just one particular operator. 
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